home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_02 / pollard.cpp < prev    next >
Text File  |  1989-04-17  |  3KB  |  99 lines
  1. /*
  2.  *   Program to factor big numbers using Pollards (p-1) method.
  3.  *   Works when for some prime divisor p of n, p-1 has only
  4.  *   small factors.
  5.  *   See "Speeding the Pollard and Elliptic Curve Methods"
  6.  *   by Peter Montgomery, Math. Comp. Vol. 48 Jan. 1987 pp243-264
  7.  */
  8.  
  9. #include <stream.hpp>
  10. #include "number.hpp"
  11.  
  12. #define LIMIT1 1000        /* must be int, and > MULT/2 */
  13. #define LIMIT2 30000L      /* may be long */
  14. #define MULT   210          /* must be int, product of small primes 2.3.. */
  15.  
  16. miracl precision=50;
  17.  
  18. static GFn bu[1+MULT/2];
  19. static bool cp[1+MULT/2];
  20.  
  21. main()
  22. {  /*  factoring program using Pollards (p-1) method */
  23.     int phase,m,iv,pos,btch;
  24.     int *primes;
  25.     long i,p,pa;
  26.     GFn b=2;
  27.     GFn q,n,t,bw,bvw,bd,bp;
  28.     primes=gprime(LIMIT1);
  29.     for (m=1;m<=MULT/2;m+=2)
  30.         if (igcd(MULT,m)==1) cp[m]=TRUE;
  31.         else                 cp[m]=FALSE;
  32.     cout << "input number to be factored\n";
  33.     cin >> n;
  34.     modulo(n);                    /* do all arithmetic mod n */
  35.     phase=1;
  36.     p=0;
  37.     btch=50;
  38.     i=0;
  39.     forever
  40.     { /* main loop */
  41.         if (phase==1)
  42.         { /* looking for all factors of (p-1) < LIMIT1 */
  43.             p=primes[i];
  44.             if (primes[i+1]==0)
  45.             { /* now change gear */
  46.                 phase=2;
  47.                 bw=pow(b,8);
  48.                 t=1;
  49.                 bp=bu[1]=b;
  50.                 for (m=3;m<=MULT/2;m+=2)
  51.                 { 
  52.                     t*=bw;
  53.                     bp*=t;
  54.                     if (cp[m]) bu[m]=bp;
  55.                 }
  56.                 t=pow(b,MULT);
  57.                 t=pow(t,MULT);
  58.                 bd=t*t;
  59.                 iv=p/MULT;
  60.                 if (p%MULT>MULT/2) iv++,p=2*(long)iv*MULT-p;
  61.                 bw=pow(t,(2*iv-1));
  62.                 bvw=pow(t,iv);
  63.                 bvw=pow(bvw,iv);
  64.                 q=bvw-bu[p%MULT];
  65.                 btch*=10;
  66.                 i++;
  67.                 continue;
  68.             }
  69.             pa=p;
  70.             while ((LIMIT1/p) > pa) pa*=p;
  71.             b=pow(b,(int)pa);
  72.             q=b-1;
  73.         }
  74.         else
  75.         { /* phase 2 - looking for one last prime factor of (p-1) < LIMIT2 */
  76.             p+=2;
  77.             pos=p%MULT;
  78.             if (pos>MULT/2)
  79.             { 
  80.                 iv++;
  81.                 p=(long)iv*MULT+1;
  82.                 pos=1;
  83.                 bw*=bd;
  84.                 bvw*=bw;
  85.             }
  86.             if (!cp[pos]) continue;
  87.             q*=(bvw-bu[pos]);
  88.         }
  89.         if (i++%btch==0)
  90.         { /* try for a solution */
  91.             t=gcd(q,n);
  92.             if (t==1) continue;
  93.             cout << "\nfactor is\n" << t << "\n";
  94.             exit(0);
  95.         }
  96.     } 
  97. }
  98.  
  99.